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ABSTRACT 

We present a comprehensive methodology for the simulation of astronomical images from optical 
survey telescopes. We use a photon Monte Carlo approach to construct images by sampling photons 
from models of astronomical source populations, and then simulating those photons through the 
system as they interact with the atmosphere, telescope, and camera. We demonstrate that all physical 
effects for optical light that determine the shapes, locations, and brightnesses of individual stars 
and galaxies can be accurately represented in this formalism. By using large scale grid computing, 
modern processors, and an efficient implementation that can produce 400,000 photons/second, we 
demonstrate that even very large optical surveys can be now be simulated. We demonstrate that we 
are able to: 1) construct kilometer scale phase screens necessary for wide-field telescopes, 2) reproduce 
atmospheric point-spread-function moments using a fast novel hybrid geometric/Fourier technique for 
non-diffraction limited telescopes, 3) accurately reproduce the expected spot diagrams for complex 
aspheric optical designs, and 4) recover system effective area predicted from analytic photometry 
integrals. This new code, the photon simulator ( PhoSivn ), is publicly available. We have implemented 
the Large Synoptic Survey Telescope (LSST) design, and it can be extended to other telescopes. 

We expect that because of the comprehensive physics implemented in PhoSim, it will be used by 
the community to plan future observations, interpret detailed existing observations, and quantify 
systematics related to various astronomical measurements. Future development and validation by 
comparisons with real data will continue to improve the fidelity and usability of the code. 

Keywords: atmospheric effects- telescopes- instrumentation: detectors- surveys- Galaxy: structure- 
cosmology: observations 


1. INTRODUCTION 

1.1. Astrophysical Instrument Simulations 

Modern survey astronomy is leading to extremely large 
samples of source populations amenable to subtle statis¬ 
tical analyses. The precision astrophysical investigations 
that are thereby enabled require high fidelity simulations 
to properly interpret the observations. In addition to 
detailed simulations of astrophysical environments (e.g. 


Fryxell et al.pOOO Springel|2005 Stone et al.|2008| ), it is 

increasingly important to build nigh fidelity end-to-end 

simulations of the instrumentation as well. This can be 
important for design and optimization of the telescope 
and camera, for planning future astronomical observa¬ 
tions, and for building sophisticated analysis software 
prior to data taking. In addition, simulated data can be 
processed along with the actual observed data through 
the same analysis pipelines, so that the flaws, biases, and 
efficiencies of the analyses can be determined accurately. 
Such an approach is now an indispensable part of mod¬ 
ern p article physics (Agostinelli et al. 2003; All ison et al. 
2006), high energy astrophysics (|Peterson, Jernigan, fe ! 
ahn||2004t |Peterson, Marshall, fe Andersson||2007 ' An- 


dersson, Peterson, fe Madejski|[20071 |Davis et al.]12012[ 

Ackermann et aT|2012), optical astronomy using adap- 


five optics (Lane, Glindemann & Dainty 1992 

Ellerbroek 


2002 Le Louarn 2002 Britton 2004 Jolissaint 

2010), ana 

has recently become more common in optical survey as- 


tronomyJ|Bertin 

2009 

Dobke et al. 

2010 

Mandelbaum 

et al. 2012). In this wor 

k, we outline a method lor high h- 


deiity optical astronomical image simulation appropriate 
for survey telescopes. 


1.2. Goals for the Simulation 

The goal of this work is to produce high fidelity sim¬ 
ulated optical astronomical images. The level of detail 
necessary to achieve high fidelity can be precisely de¬ 
fined by considering the measurable properties of the 
images that we are interested in reproducing. We can 
divide those measurable image properties into two cat¬ 
egories: 1) primary image properties and 2) secondary 
image properties. The primary image properties are: the 
point-spread-function (PSF) full width at half maximum 
(FWHM), the photometric zeropoint, the plate scale (or 
astrometric scale), and the background intensity level. 
These quantities are predicted from physical effects in 
the instrument/atmosphere system. If we are only inter¬ 
ested in simulated images that correctly reproduce those 
four primary properties, it is possible to use an empirical 
parametric approach. For example, the photometric ze¬ 
ropoint and the object’s flux within a band can be used to 
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calculate an object’s intensity. The PSF FWHM and the 
object’s spatial distribution yield its observed morphol¬ 
ogy. The astrometric scale determines the location in the 
image, and the background provides the extra uniform 
intensity level. In many cases, such a simulator would be 
sufficient to predict what a basic image would look like 
from an astronomical telescope. However, a simplified 
simulation tool that only matches these four properties 
is insufficient for detailed applications. 

It is significantly more difficult to predict more com¬ 
plex secondary image properties using parametric mod¬ 
els. Examples include: the PSF radial profile, the PSF 
size wavelength dependence, the spatial variation of the 
PSF size, the PSF shape (notably, ellipticity), the PSF 
shape wavelength dependence, the PSF shape spatial 
decorrelation, the PSF shape spatial variation, the differ¬ 
ential astrometric shift non-linearity with field angle, the 
differential astrometric non-linearity with wavelength, 
the differential astrometric decorrelation with angle, the 
differential astrometric decorrelation variation, the pho¬ 
tometric chromaticity, the photometric variation in time, 
the photometric variation with field angle, the back¬ 
ground variation in time, the background spatial depen¬ 
dence, and the background wavelength dependence. Un¬ 
derstanding such properties is essential for subtle statisti¬ 
cal analyses associated with galaxy shape measurements, 
stellar astrometry, and precision photometry. Some of 
the astrophysical investigations that require such anal¬ 
yses include: photometry of superno vae for measuring 
the expansion rate of the Universe (Riess et al. 1998 
Perlmutter et al. 1999), using stellar proper motions ana 


parallax measurements to map the structure and kine¬ 


ma tics of the Milky Way (Hoeg et al. 2000 Monet et 


al. 

2003 

1, 

measurement of 

large 

scale galaxy structure 

sta 

tistics 

Eisenstein et al. 

2005 

), and studies of weak 


gravitational lensmg for bo 


h dark matter and dark en¬ 


ergy investigations (Tyson, Wenk, ~fc Valdesf l990). 

In this paper, we describe a simulation tool that takes 
explicit account of all the detailed atmosphere, telescope, 
and camera physical effects that determine the propaga¬ 
tion of light to the focal plane. We show that the relevant 
physics can be most efficiently encoded in terms of pho¬ 
ton manipulations. The goal of this work is to be able to 
predict all of the primary and secondary image properties 
from physical models of the atmosphere and instrument 
within the context of a photon Monte Carlo approach. 
We use very few approximations, and the overall simu¬ 
lation accuracy has been carefully controlled to ensure 
that calculational errors are small in comparison to the 
statistical uncertainties associated with measurements of 
typical sources in the survey. 


1.3. Simulation Regime for Optical Surveys 

In this work, we are interested in simulating images 
containing large numbers of astronomical objects (typi¬ 
cally, millions) in the visible and near infrared wave-band 
(^300 to ^1200 nm) that would be obtained in an opti¬ 
cal survey. A complete physical simulation with no ap¬ 
proximations would take account of the detailed material 
properties of every element of the telescope and camera 
and their response to all external forces and tempera¬ 
ture variations, a complete hydrodynamic treatment of 
the atmosphere, and simulation of the full quantum me¬ 
chanical nature of light. However, such a rigorous ap¬ 


proach is unwieldy, and unnecessary, since a number of 
key simplifying approximations can be introduced for op¬ 
tical survey telescopes which are non-diffraction limited. 
A telescope is non-diffraction limited if 


max 


A a 
ro’F 


>> 


A 

D 


where A is the wavelength, D is the diameter of the tele¬ 
scope, F is the focal length of the telescope, ro is the 
Fried parameter, and a is the characteristic size of aber¬ 
rations introduced by the telescope and camera. The 
Fried parameter describes the characteristic beam diam- 
eter w he re t he optical phase rms value is close to 1 radian 
(Fried|l965). Thus, the first ratio is associated with dis¬ 
tortions introduced by optical turbulence (usually named 
seeing), while the second accounts for imperfections in 
the optical system. 

It is convenient to rewrite this condition in the form: 


max 


D a 
n>’ A / 


» 1 


where / is the focal ratio of the system, F/D. For most 
ground-based optical telescopes the value of ro lies be¬ 
tween 5 and 30 cm, depending on atmospheric condi¬ 
tions. For modern survey telescopes like LSST (8.4 m 
diameter primary), the first ratio is much larger than 
unitj0 In addition, to achieve wide fields of view over 
a range of wavelengths, practical issues cause these tele¬ 
scopes to be somewhat aberrated, so that the second ra¬ 
tio is generally greater than unity as well. If some of the 
physical effects (e.g. either large aberrations from atmo¬ 
spheric turbulence or large telescope optics aberrations) 
cause the system to be non diffraction-limited, we can 
consider standard ray optics techniques for those parts 
of the photon simulation. This does not preclude the 
use of quantum mechanical calculation or wave-like in¬ 
terference computational techniques at the appropriate 
physical places as we discuss in §2.1. 


1.4. LSST 

The initial implementation of the data describing 
a telescope and site in PhoSim is LSST. LSST is a 
wide-field ground-based telescope with extremely large 
etendue (the product of the effective area and the field 
of view), 320 m 2 deg 2 , roughly a factor ten higher than 
all previous facilities. The large etendue will enable an 
unprecedented survey of the optical sky. It also results 
in a very high anticipated data rate of 15 Terabytes of 
images per night. LSST is an extreme case of the optical 
survey regime discussed above. It is also the most chal¬ 
lenging telescope to simulate computationally because of 
the large quantity of data. We have therefore chosen to 
construct the simulator to simulate LSST and focus on 
that application primarily in this paper. However, we 
have written the simulation code to keep the LSST de¬ 
sign data separate from the physics algorithms, so other 
telescopes can be simulated with this code as well. 

1 6 

The Fried parameter depends on wavelength as As, so 

throughout this work, values are referenced to 500 nm as this is 
the standard convention in the literature. 
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LSST will be constructed on the El Pehon peak at the 
Cerro Pachon ridge in Chile at 2660 m above sea level. 
The site is expected to deliver a median seeing of 0.67 
arcseconds (with 0.44” as the 25tlr percentile and 0.81” 
as the 75th percentile). LSST is an 8.4 m optical tele¬ 
scope with a three mirror modified Paul Baker design 
to allow for a large unaberrated field of view. It has an 
effective f/# of 1.23. It has three correcting lenses and 
complement of filters as part of a 3.2 gigapixel camera. 
The camera has 189 individual 4k x 4k CCDs. The 6 fil¬ 
ters (u,g,r,i,z,y) cover the wavelength region from 300 to 
1100 nm. The intrinsic aberrations of the optical system, 
the residual aberrations after compensation of the active 
optics control system using curvature sensors, and the 
charge diffusion in the devices is expected to contribute 


about 0.4 arcseconds to the PSF (|Ivezic et al. 2008 LSST 
Reference Design||2012 ). 


1.5. Scope and Interfaces 

The simulation work described here is embodied in a 
publically available code called the Photon Simulator or 
PhoSim. The code and the documentation are located 
at https: / / www. bitbucket. org/phosim / phosim_release. 
PhoSim is a stand-alone code that requires catalogs of 
astropliysical objects (positions, fluxes, shape, spectra, 
etc.) and operational parameters (pointing information, 
instrument configuration, etc.) as input. The output 
of PhoSim is a stream of images. The user can create 
the inputs to either match a real observation or to 
create a hypothetical observation. The combination of 
the operational parameters and astropliysical catalog 
data is called an instance catalog and is described in §3. 
Internally, PhoSim uses the physics of the atmosphere, 
telescope, and camera to predict high fidelity images 
based on the input and is described in the following 
section. 


2. PHYSICS DESCRIPTION AND ALGORITHMS 

In the following, we describe the complete set of at¬ 
mosphere/instrument physics that is encapsulated in the 
PhoSim code. We describe how the physical interactions 
can be computed in terms of photons manipulations. 


electromagnetic wave interference calculation at interface 
boundaries. Thus, ray optics is the methodology we use 
to represent the photon’s interaction with physical el¬ 
ements where ray optics is appropriate, whereas more 
complex wave optics (or quantum mechanical physics) 
can describe the photon in other areas. 

In ray optics, the photon state can be represented by 
a set of eight numbers (ignoring polarizatiorrl) : a vec¬ 
tor describing the photon’s current position (x), a unit- 
vector describing the angle of propagation (n), a time 
stamp when the photon arrived at Earth (t), and a wave¬ 
length for the photon (A). Throughout the simulation of 
the physics, 3 possible alterations of the photon’s vari¬ 
ables may occur: 1) the photon’s trajectory may change 
and therefore the unit-vector changes its value, 2) the 
photon’s position may change as it propagates along its 
path, and 3) the photon may be removed (most likely 
because it scattered in an uninteresting direction). To a 
good approximation, the photon’s wavelength does not 
change and the time stamp of the photon can be re¬ 
garded a constant, since a photon moves from the top 
of the atmosphere to the detector in less than 30 ps 
and none of the physical models have rapid changes on 
that time-scale. By the end of the calculation, we may 
have converted the photon into an electron in the detec¬ 
tor. The electron similarly can be described by the same 
formalism, except that the magnitude of the velocity is 
the relevant quantity instead of the wavelength. We are 
therefore using a photon Monte Carlo method to simul¬ 
taneously manage the change in trajectory as well as the 
change in intensity or throughput. 

2.2. Sky Photon Sampling 

The first step of the simulation is to create the photons 
from astrophysical sources. This involves populating the 
quantities describing a photon. The propagation direc¬ 
tion for a photon coming from a source at a position 
( ai,5i ) is determined by a unit vector, £,(/, and z that 
is calculated from the current bore-sight of the telescope 
(a, S) according to 

z = (cos a cos 6 , sin a cos 6, sin 6) 


2.1. Basic Monte Carlo Methodology 


Given the simulation regime (§1-3) and the simula¬ 
tion goals (§1.2) it is ideal to adopt a Monte Carlo ap¬ 
proach. The basic Monte Carlo approach has been an 
essential part of computational physics for some tim e 


( Ulam et al.|1947j [Metropolis fc Ulam|1949 Ulam|1 950[). 
In our case, the Monte Carlo is implemented by l 


implemented by follow- 


ing individual astrophysical photons. Thus the integral 
over the time-dependent, angle-dependent, and energy- 
dependent incident radiation field is accomplished by 
simply following individual photons throughout the at¬ 
mosphere/instrument system. This is particularly im¬ 
portant because a large fraction of the physical effects 
have wavelength and angle-dependent properties. 

In certain contexts, the photon’s propagation can be 
represented using ray optics. However, as we show in this 
formalism, this does not preclude us from using a diffrac¬ 
tion calculation for perturbing elements with a wavefront 
shift smaller than the photon wavelength, a quantum me¬ 
chanical calculation for photon-atom interactions, or an 


x = (cos ( a + 7r/2), sin (a + 7r/2), 0) 


y = z x x 

where x represents a vector cross product. The coordi¬ 
nate system is then in the telescope’s frame where z is 
the optical axis. We calculate the angle of propagation 
by first calculating a source vector, s 

s = (cos OLi cos Si , sin a* cos Si , sin Si ) 
and convert the vector into tangent coordinates by 

s-x „ s-y „ r — — 

% = — ; n 9 = — ; n z = -J1 - n 2 - n 2 

s ■ z s ■ z V y 

where • indicates a dot product. The source’s position is 
assumed to have astronomical aberration included. The 

2 Polarization is relevant in interference calculations that result 
in fringing. 











4 


Peterson et al. 


time stamp, t of the photon is chosen in a Poisson manner 
by considering the total exposure time, t e and a uniform 
random number, u 

t = to + t e u 

where we include a time offset, foi if we are simulating a 
sequence of exposures. This will affect other quantities 
in the case of transient or moving objects. 

We determine the photon’s initial position, x, by first 
sampling the annular pupil of the primary mirror by 
choosing the photon’s position in polar coordinates, 



where u and v are uniform random numbers, r Q is the 
outer radius and r, is the inner radius. The outer and 
inner radius can be chosen to be slightly larger than the 
actual physical aperture to allow for photons that scatter 
into the aperture. Then, we calculate x = r cos (f> and y = 
r sin (f >, and determine the z position by using the mirror 
surface function, z(x,y). Then, we move the photon to 
the top of the atmosphere x x + nl. I is calculated 
using — h 'l~ z where ho is the top of the atmosphere (100 
km). 

We choose the wavelength of the photon by sam¬ 
pling from a spectral energy distribution describing the 
source. Spectral energy distributions (SEDs) are nor¬ 
mally quoted in terms of flux units (ergs cm -2 s _1 A ) 
at given wavelength interval in the spectrum. We can 
convert this into a relative probability of a photon ap¬ 
pearing in a given wavelength interval by dividing by 
the average photon energy and multiplying by the wave¬ 
length interval. Then we sample from this distribution 
by converting into a cumulative distribution and draw¬ 
ing a uniform random number. We blur the wavelength 
by the bin width. By convention, we use the part of the 
SED with wavelength less than 1200 nm (which would 
be appropriate for LSST). We also draw the photon in 
the rest frame of the object, and then redshift the pho¬ 
ton by multiplying the wavelength by 1+z. This allows 
us to re-use the SEDs for different galaxies at different 
redshifts, and still obtain a diversity of galaxy colors. 


ar/ro and br/ro where a and b are model inputs. Fi¬ 
nally, we rotate the major and minor axes by a rotation 
angle, <f >. In practice, we have found it useful to represent 
galaxies as a pair of ellipsoidal Sersic models where one 
represents the bulge and one represents the disk. Each 
component has a different Sersic index, their centers are 
possibly offset, and they can have different major and 
minor axes. A pair of ellipsoidal Sersic models results in 
a 12 parameter model with reasonable fidelity. 

For more accurate galaxy morphology simulations, an¬ 
other common model we include is to simply input truth 
images of the source before it would have been observed 
through the effects of the atmosphere and telescope. We 
use this by using this truth image as the 2 dimensional 
probability distribution of finding a photon in a given 
angular pixel. We then select photons from this distribu¬ 
tion. We also allow for an arbitrary rotation and scaling 
that places the source at different distances and orienta¬ 
tions. The truth images can either come from other tele¬ 
scopes (particularly those taken with a night with better 
seeing or a space-based telescopes), or be generated by 
another simulation. 

A simple perturbation to the spatial models is to in¬ 
clude the effect of weak gravitational lensing for distant 
galaxies. This is accomplished by taking the photon’s rel¬ 
ative position to the source’s center ( Sx , 5y) and apply 
the matrix 

(l-li “« 72 \ 

\ 72 1 + 7i ~ K ) 

where 71 , 72 , and k are the usual weak lensing shear and 
convergence parameters. The user’s cosmological simu¬ 
lation can estimate these parameters from the individual 
foreground dark matter distribution for each galaxy in¬ 
dividually. We handle moving objects (e.g. Solar System 
objects or satellites) by simply perturbing the initial pho¬ 
ton direction by a proper motion vector (in units of angle 
per time) and multiplying by the photon’s relative time 
stamp. This results in streaks when the full simulation 
is completed. Both of these effects can also be handled 
by a distorted truth image, but are added options since 
they can be applied to objects when expressed as catalog 
entries. 


2.2.1. Non-Point Source Models 

Sources that are extended on the sky are treated by 
first choosing the direction, n, using the procedure above 
for the nominal center (a*, S t ) of the source emission, but 
then perturbing the unit vector by a structural model for 
the spatial emission. We can use a paraxial ray approx¬ 
imation to perturb the values of h x and n y and then 
renormalize the unit vector. 

The most useful spatial model we have constructed is 
an ellipsoidal Sersic distri bution. The intensity of light 
of the Sersic distribution (Sersic] 1963) is given by 


J(r) = Joe b "( r o) 

where ro is the scale radius, n is the Sersic index, and 
b n is a normalization constant for each n. We then draw 
a value of r from this probability distribution, and then 
choose the distance along the major and minor axes as 


2.2.2. Diffuse Emission Simulation 

There are a large number of photons that contribute 
to the optical sky background. We model these pho¬ 
tons from 3 sources: airglow from the dark sky, reflected 
moonlight, and additional emission near twilight. We 
also simulate light from a dome screen through the same 
mechanism. For each of these forms of diffuse light, we 
model them by a collection of sources uniformly spaced 
on the celestial sphere. This is necessary to allow for 
complex spatial gradients. We space the sources by 15 
arcseconds, and draw photons from a two dimensional 
spatial gaussian having a standard deviation of 15 arc- 
secondf£] This results in a statistically uniform illumina¬ 
tion pattern, but limits us to not have a spatial variation 
of the illumination patterns on scales smaller than 15 

3 For arbitrarily high statistics this approach will have a periodic 
wave artifact, however, for practical background levels the spacing 
of the sources is sufficient to make this unobservable. 
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arcseconds. For each diffuse model, we therefore need 
to predict the intensity of light across the field on 15 
arcseconds scale. 

For the simulation of airglow emission , we use a spec - 
tral energy distribution taken from Patat et al. (2006). 
The overall intensity for a given exposure is simulated 
as having a r band magnitude from a gaussian distri¬ 
bution with mean of 22.08 and a standard deviation of 
0.9. We then predict a relative variation across th e fiel d 
using spatia l power spectrum measurements of [Ad ams| 


& Skrutskie (19961 where the variation was proportional 
to fc _oT "and normalized to be 3% at 1 degree. We also 
in crease the emission in pr oportion to the zenith angle 
(Adams & Skrutskie 1996). The mean sky background 
a nd its variatio n are consisten t with the measurements 
of Krisciunas & Shaefer (1991). 

The moon’s intrinsic brightness as a fun ction of its 
phase and alti tude follows the calculation oflKrisciunas fe 


Shaefer (1991). We use an empirical lunar spectrum for 
its SED. We then need to predict the bright ness where 
the telescope i s pointing. Here we use the Krisciunas fe| 


Shaefer (19911 formula that has ter ms for the Rayleigh 
and Mie s cattering of the moonlight. |Krisciunas fe ShatT 


fer (1991) only calculated the lunar brightness for one 


band so we simply scale the Rayleigh term by inverse 
wavelength to the fourth power. Conversely, Mie scat¬ 
tering is approximately wavelength independent. The 
sky brightness is increased near twilight according to the 


sun ’s alti tude using a color-dependent model of Patat et 


al. (2006). 


2 . 2 . 3 . Number of Photons 

The total number of photons for a particular source is 
calculated by considering the AB magnitude at a par¬ 
ticular wavelength. We have found it most convenient 
to normalize the spectrum at a wavelength of 500 nm 
divided by 1 + z of the source. We then convert the 
flux nearest that wavelength to units of ergs cm -2 s -1 
Hz -1 . Then, we convert the SED to a relative fraction 
of photons in each bin. Using the probability of finding 
a photon in the bin near the reference wavelength, 500 
nm/(l + z), one can then calculate the total number of 
photons per sq. centimeter per second from that source 
at all wavelengths. This convention is most useful us¬ 
ing a un-redshifted spectrum, since the same SED can 
be used for multiple sources at different redshifts. Con¬ 
versely, this does not preclude a user redshifting a SED 
before input and then setting the redshift to zero. The 
total number of photons (without including any efficiency 
losses) is calculated by multiplying by the aperture area 
and exposure time. We then modulate that expected 
photon count rate using a Poisson distribution. 


2 . 2 . 4 . Dust Models 

A final step of the sky simulation is to remove some 
fraction of the photons due to dust absorption. We use 


two du st extinction cu rves: the Cardelli, Clayton, Mathis, 
(1989) model and the Calzetti et al. (2000j) model. Each 
model has two parameters: Ay and K v . The models are 
stored in a variety of grid lookup tables and the dust is 
applied by first calculating the optical depth, r, using 
the photon’s wavelength. Then we destroy photons with 
probability e~ T . Performing dust extinction through a 


Monte Carlo approach conveniently also avoids construc¬ 
tion of a unique SED for every single source. For extra- 
galactic sources, we have the option of applying the dust 
absorption both before and after the redshift of the pho¬ 
ton, representing absorption in the galaxy as well as in 
the Milky Way. 

2 . 3 . Atmosphere Simulation 

After creation of the photons from the astrophysical 
and sky emission, we then propagate the photons through 
an atmosphere simulation. We first describe the struc¬ 
ture of the atmospheric components, and then describe 
how the photon interacts with those components. 


2 . 3 . 1 . Atmosphere Structure 

We model the atmospheric structure by a series of 
plane-parallel layers. Each layer covers the absorption 
and turbulence blurring between two altitudes. The 
propagation of a ray from one layer to the next is 
straightforward using plane-parallel layers. The photon’s 
position can be updated using x —> x+hl where the scalar 
distance, l is calculated from 


l = 


hi — h 


i +1 


where hi is the altitude of the interface between the two 
layers and h z is the same as in §2.2. 

Atmospheric turbulence is known to significantly de¬ 
grade overall image quality of astronomical telescopes. 
The energy power spectrum of atmospheric turbulence 
is known to have a power spectrum that approximately 
follows the scale-free Kolmogorov spectrum (A -11 / 3 ), so 
the largest powers are on the largest scales ([Kolmogorov 
19411. This spectrum extends down to the viscous limit 
(a few mm) and up to the scale where the turbulence is 
driven. The complete spectrum is often parameterized 
by an oute r scale , Lq, an d a f lat power spectrum above 
that point (von Karnian] 1948]), 


(k 2 + L^)-" le 

The justification of using this spectrum to describe the 
real at mosp h eric turbul ence has been est ablished in de¬ 
tail (Tokovinin, Sarazin, & Smette 2007). The use of a 
von Karman spectrum rather than a Kolmogor ov spec¬ 
trum has several observational consequences (Martinez 


et al. 20101. Other more complex power s pectra are 
sometimes considered (Hill & Clifford 1978) and non- 
Kolmogorov effects are an active area of turbulence re¬ 
search, but astronomical image quality is generally not 
particularly sensitive to these details. 

The temperature variation produced by this airflow 
turbulence pattern leads to index of refr action va r iation s 
that affect light propagation (see e.g. Roddier 1981). 
The deviation of index of refraction from unity (n — 1) is 
approximately linearly proportional to the pressure and 
inversely to the temperature, and the effect of pressure 
fluctuation is negligible compared to the temperature 
variation. So therefore the same turbulent eddies pro¬ 
duce index of refraction variations with the same spec¬ 
trum as the temperature variation. The entire effect of 
the turbulence on light propagation is then represented 
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by a phase shift as a function of spatial position (i.e., a 
phase screen). 

Furthermore, the turbulence pattern itself for a par¬ 
ticular layer is known to not change significantly during 
the time it takes to cross the aperture of the telescope 
(a fraction of a second). To a good approximation, the 
pattern tends to drift with some wind velocity vector 


ing (Taylor hypothesis: 

Taylor 1938 

Favre, Gaviglio, & 

Dumas 1952 Foyneer et al.|2009). In addition, the most 


atively narrow interfaces due to either differential shear¬ 
ing of the atmosphere or vertical instabilities. Hence, it 
is very common to represent the turbulence as a series of 
frozen plane-parallel two-dimensional screens where the 
three-dimensional structure for a given slab of the atmo¬ 
sphere has already been collapsed into two dimensions. 

We therefore propagate light through a series frozen- 
phase screens drifted by wind velocity vectors. This 
simulation approach is standard in the adaptive op- 


mann & Dainty 

1992, 

Ellerbroek||2002, 

Le LouarnP002, 

Britton 2004, Joiissaint 2010|. Howevt 

3r, we made two 

problem that are 

novel innovations to suit our particular 


not standard techniques. First, instead of propagating 
light through by computing the full diffraction integral, 
we used a novel geometric raytracing approximation for 
the low-frequency part of the phase screen described in 
§2.3.2. Second, we constructed the phase screens on 
four different scales and repeated the three smaller phase 
screens on the larger s cales. The tiling s cheme has been 


previously explored by Vorontsov et al. (20081 for a sin¬ 
gle scale, and here we simply extend the method to four 
scales. This general approach is necessary to capture the 
large scales involved in the simulation of a wide-field tele¬ 
scope like LSST. In particular, the screen size, L has to 
be large enough to not repeat during an exposure time 
(L > ^ w i nc ^exp) and large enough to not repeat when 
off-axis sources are considered (L > h9 , where h is the 
height of the highest layer). We therefore use a linear 
superposition of four 1024 by 1024 pixel screens having a 
pixel size of 1, 8, 64 and 512 cm. Every pixel in the 5 km 
x 5 km screen is therefore unique. Numerical artifacts 
result from certain alignments of patterns along the fun¬ 
damental axes of the screens with the wind direction (0, 
f, arctanf, etc.) since the patterns may be repeat 
and a small part of the screen gets used. However, the 
actual atmosphere has about a 5 degree per minute drift 
of the wind direction (Mahrt 2010), which makes these 
artifacts not occur in the complete simulation when this 
effect is turned on. This is particularly important in the 
detailed simulations of wavefront images. 

We drift the turbulence screens according to a wind 
model where we predict a wind vector as a function of 
height for a particular observation. In order to deter¬ 
mine where a photon hit a screen at a given layer, we 
first calculate the x and y position. We then calculate 
the pixel in the appropriate screen given two components 
of the wind vector for each screen. The arrival time of the 
photon then dictates exactly which pixel is used. Using 
the NOAA NCEP/NCAR Reanalysis Monthly Database 
for a particular location, we fit the historical data to a 
Weibull distribution (k=2) for the wind velocity. Both 


the wind direction and magnitude have a seasonal distri¬ 
bution. We drift the wind dir ection during an exposure 
according to measurements of Mahrt (2010). 

In order to construct the turbulence screens, we require 
the single parameter, the outer scale, of the |von Kar~7| 
man] (1948) model for each screen. Within the context of 


this model, it is often argued that this varies as a func¬ 
tion of height, has a large time variation, and likely has 


a log-normal distribution (e.g. Beland fe Brown 1988 


Coulman et al. 1988] Abahamid et al.| 2001 [) . A log¬ 


normal distribution is a good match to the LSST site 
where the mean (35.6 m) an d m edian (26.7 m) param¬ 
eters were measured by Boccas (]2004|). The possible 
altitude-dependence of the outer scale is still an active 
area of research (e.g Abahamid et al. 2004), so we do not 
have an altitude-dependence to the outer scale. 

To predict the relative turbulence intensity as a func¬ 
tion of altitude, we construct an interface based on the 


model of Tokovinin & Travouillon (2006) for the Cerro 
Pachon site. However, other sites can clearly be repre¬ 
sented by different set of parameters as is evident in the 
general nature of the parameterization below. In this 
model, they use a 7 layer model (a ground layer and 6 
layers at 0.5, 1, 2, 4, 8, and 16 km). With a series of mea¬ 
surement they then quantified the 25th, 50th, and 75th 
percentile of the turbulence integral, J = f C%(h)dh, 
where Cfj is the refractive index structure function and 
h is the altitude. We then use these percentiles and rep¬ 
resent the complete distributi on of J at each altitude 
as a lognormal distribution. Tokovinin & Travouillon 
(2006) found that the turbulence intensity of the ground 
layer was largely independent of the other layers. The 
free atmosphere layers, however, were highly correlated 
as seen by comparing Figure 5 and Table 3 in Tokovinin 
& Travouillon| (2006). Therefore, we use a single ran¬ 
dom number for the free atmosphere and a single ran¬ 
dom number for the ground layer. This then tends to 
produce poor seeing when either the free atmosphere or 
the ground layer is particularly turbulent. We can use 
this model in two different ways. First, if a random atmo¬ 
sphere is desired, we simply draw turbulence intensities 
from this model. Second, to obtain a particular seeing, 
we run this model hundreds of time to produce the de¬ 
sired seeing. In either case, the final turbulence intensity 
values are used to normalize each screen. 

In addition to the screens of turbulence, we track four 
components that represent the bulk air located between 
the screens: the density of the atmosphere vs. height, 
the density of water vapor vs. height, the density of 
molecular ox ygen vs. height , an d the density of ozo ne 


vs. height (Sander et al. 2006 Thomas & Stramnes 


19991. These profiles are the n used~t o calcul ate the col¬ 


umn density between layers (Bodhane et al. 1999 


Liou 


2002) We calculate the column density relative to tne ap¬ 
propriate altitude. We also modulate the overall column 
density from one simulation to the next by a lognormal 
distribution with width equal to 0.18, 0.20, 0.002, and 
0.01 for water, ozone, molecular oxygen, and the overall 
density, respectively. This represents the natural time- 
dependent variation. When the telescope is not pointed 
at zenith, we increase the column density by the exact 
zenith-dependent airmass. We recalculate this for every 
photon separately since the variation of airmass across 
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the field can be significant for wide fields. We then use 
these column densities to calculate the optical depths, 
which is described in §2.4.2. The local column density is 

also perturbed by a factor of p x /./Vi of a screen 

1 J 1 A/ layers 


with an arbitrary variation pattern at each layer. In 
this way, the opacity will vary slightly from exposure to 
the next, and it will vary across the field in a complex 
way with amplitude, p. However, the parameterization 
of reasonable variation patterns is still an active area of 
research. We also have clouds screens to represent the 
relative opacity of cloud absorption at various heights. 
The cloud screens have an exponential structure function 
with a angular coherence scale of 2 de grees to be consis¬ 
tent with the measurements of Ivezic et al. (2007) based 
on SDSS wide field images. We assumed that the struc¬ 
ture is isotropic, but clearly have non-isotropic structure 
that could be related to the wind direction. We use two 
layers of clouds each with a different wind velocity to 
create a partly realistic complexity that would at least 
mimic the difficulties for photometric calibration. The 
photon interactions with all the components and prop¬ 
erties of the atmosphere are described in the following 
sections. 


2.3.2. Photon Interactions with the Turbulent Atmosphere 


An important component of the atmosphere simulation 
is the propagation of the photons through the turbulent 
screens. The theory of how light propagates t hrou gh a 


turbulent medium has been extensively studied ( 

Tatarski 

1961 

Fried 1965, Clifford 1978, Roddier 1981 

Schmidt 

2010 

difir; 

1 , and it is common to use Hugyens-Fresnei scalar 
iction theory, which represents light as a monochro- 


rnatic wave having a scalar amplitude at all spatial po¬ 
sitions. The so-called geometric optics limit , where non¬ 
linear effects can be ignored, is valid when the charac¬ 
teristic height of turbulent layers in the atmosphere is 
less than the square of the characteristic turbulent cell 
size divided by the wavelength in which the observation 
is performed 


, 12 

hi j 

For typical atmospheric conditions at good sites, h ss 10 
km, and 1 ss 10 cm, which means that this condition is 
essentially satisfied across the entire optical band: A = 
0.3 to 1 tun. In this limit, we can write th e PSF in term s 
of the Fraunhofer integral (see e.g. [Goodman |[20Of]). 
Here, the instantaneous PSF is given by the square of 
the Fourier transform of the electric wave amplitud^] 
over the pupil plane 


PSF{n) 


1 

A 


fpfeikn-r^lr) 


where k = 27r/A, n is the unit vector pointing to that 
particular point on the focal plane, and A is the tele¬ 
scope pupil transmission. For the usual small angular 


4 The electric wave amplitude, U(r), can be written as U(r) = 
However, for the exposures considered here the 
amplitude fluctuation is negligible and the last factor has no effect 
on the image, so only the second factor is relevant in the above 
integral. 


displacements appropriate to the PSF, h has coordinates 
( 9 X , 9 y , 1), where 9 X and 9 y are coordinates in the focal 
plane with units of angle. The phase, <fi, which appears 
in this expression is the total phase shift produced by all 
of the atmospheric layers along the line of sight due to 
variations in the index of refraction. Note that the PSF 
can also be expressed as the Fourier Transform of the 
correlation function of the exponential of the complex 
phase 


B{p) = ^ f dPfeWe-W+d 

A J a 

The turbulence in the atmosphere is a stochastic process, 
which means that both the phase correlation function 
and the PSF will depend strongly on time through the 
variation in the random structures of the various atmo¬ 
spheric layers that cross the pupil plane at any particular 
instant. However, for an exposure of finite duration, a 
larger footprint on these atmospheric layers is observed, 
which means that the effective PSF is averaged over such 
structures, and it is meaningful to estimate it via statis¬ 
tical techniques. In the adaptive optics literature, it is 
common to talk about two distinct limits: the long expo¬ 
sure limit, in which one takes a complete average over the 
turbulent structure on all spatial scales, and the short ex¬ 
posure limit, in which one takes some average over scales 
which are small compared to the diameter of the pupil, 
but ignores the effects of much larger scales. The latter 
is appropriate, because as we show below, the primary 
effect of the larger scales for short exposures is image 
displacement, not an increase in PSF width. 

If we are interested in estimating PSF anisotropy 
and decentering for moderate exposures, an intermedi¬ 
ate limit is required. Since there are no preferred axes 
in the problem, the long exposure limit, which involves 
a full statistical average over all of the turbulent struc¬ 
tures must clearly yield a circularly symmetric PSF with 
no decenter. However, there are significant variations in 
the image motion for moderate exposures, which are in 
fact the principal causes of residual anisotropies, so these 
must be modeled correctly. That means that the short 
exposure limit is also inappropriate. 

The intermediate limit can be found by separating out 
the contributions to the PSF from structures on large 
and small scales. In what follows, we provide a formal 
justification for how that works in detail. We first define 
the separate contributions to </> that are contributed by 
large and small scales 

<t>(r) =<t>>(r) + 

where 4>>(r) is the contribution from structures with 
wavenumber larger than some critical value, K cr ^, and 
ct> < ( r) is the contribution from structures with wavenum¬ 
ber less than k ci ^. We choose K cr such that structures 
with larger wavenumbers are very well sampled during 
the exposure, so that the long exposure limit is indeed 
applicable on those scales, and we can estimate their con¬ 
tribution statistically. For wavenumbers less than K cr j^, 
we effectively compute <f> < (r) directly for a simulated ex¬ 
posure by generating phase screens with appropriate res¬ 
olution, i.e. with a sampling 2tt/k ct 

Note that because the PSF is insensitive to the mean 
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phase, both </><(r) and (/>>(r) can be taken to have mean 
zero, with no loss of generality. Then, the correlation 
function, B(p), that we introduced earlier, can be written 
in the form 



Cpf? e i{4> > {r)-<t> > (r+p) e i(4> < (r)-4> < (r+P) 


As indicated, we evaluate the first complex exponential 
term in the integral by taking a statistical average of tur¬ 
bulent structures at higher wavenumbers. That clearly 
removes an explicit dependence on r, so it is possible 
to express the correlation function as the product of 
two separate correlation functions for the two separated 
wavenumber regimes 


B{p) = B > (p)B < (p) 

The PSF due to the atmosphere is the Fourier Transform 
of the correlation function. It can therefore be repre¬ 
sented as the convolution of the separate PSF contribu¬ 
tions from large and small wavenumbers 


PSF{h) = J d 2 rhPSF > (m)PSF < (m + fi) 

As indicated above, for wavenumbers larger than K cr q, 
we can average over the turbulent structure, because 
these scales are well-sampled even for moderate expo¬ 
sures. The analysis follows closely that which is typically 
invoked for the long exposure limit. 

Specifically, we assume the Kolmogorov theory of tur¬ 
bulence, which posits that energy is injected into an at¬ 
mospheric layer on large scales and cascades down to 
smaller and smaller scales via turbulent eddies, until it 
eventually can be dissipated by molecular viscosity. The 
instability is characterized by the Reynolds number. For 
typical atmospheric parameters, the Reynolds number 
remains above the critical value down to scales of order 
millimeters, while energy is injected on scales of 10 me¬ 
ters. Therefore, the structure is turbulent over a very 
broad range of wavenumbers. 


In Kolmogorov’s statistical treatment (Kolmogorov 
19411, the rate of energy per unit mass that is injected is 
equal to that which is dissipated down through the cas¬ 
cade to smaller and smaller scales. Simple dimensional 
arguments then suggest that the characteristic velocities 
obey the relation 


i r i 

v ~ e 3 L 3 

which means that the temperature fluctuations are pro¬ 
portional to is. Index of refraction variations are ba¬ 
sically proportional to the temperature fluctuations, so 
they are proportional to Li as well. Integrating over a 
layer, which is thick compared to the characteristic eddy 
size, leads to a two-dimensional phase power density pro¬ 
portional to r 3 , where r is a spatial scale on the two- 
dimensional layer. In spatial frequency space, the power 
spectral density is then proportional to k~ t . The full 
expression for this power spectral density is 

11 

$^(k) = (0.033)27t/c 2 k-t / dzCl(z) 

Jo 


Here C 2 is the so-called structure parameter. It is a 
measure of the distribution of seeing contributions as a 
function of altitude. It has units of m ~ 3 . 

We now consider B > (p). Taking the average over the 
turbulence yields 


B > (p)=\ [ #?< > 

A J A 

Expanding the complex exponential as a Taylor series, 
and recognizing that 4>>(r) is a Gaussian random vari¬ 
able with mean zero, it is clear that the odd moments 
vanish, so that B > can be written in the form 

B > (p) = e-iRAp") 

where D^p) is the phase structure function, given by 

D <p(P) =< [4>>(r + p) + <Mr)] 2 > 

Note that this can be written as 

D^p) = 2[B*(0) - B^p)] 

where 


=< <t>{r)<!>{r + p) > 

is the correlation of the phase itself. Since B^ip) is the 
Fourier transform of the phase power spectral density, so 
the phase structure function can be written as an integral 
over spatial frequency for the turbulent layer. Since we 
only want the contribution due to wavenumbers greater 
than K cr j|;, we can truncate this integral accordingly. The 
resulting expression is 


D^p) =47 t [1 - Jo(kp)]$<l,(K)KdK 

•' ,t crit 

The equations above yield a closed form expression for 
B > (p). The Fourier transform of that expression gives 
us PSF > (rh). 

For wavenumbers smaller than K cr ^, the statistical av¬ 
erage cannot be invoked if we hope to make a useful 
model of the anisotropy and decenter. In this regime, 
we have to calculate the phase screens directly. The ex¬ 
act approach would be to Fourier transform those phase 
screens, but this is computationally intensive for screens 
large enough for large field of views, especially when one 
considers the motion of the screens associated with the 
wind velocities during the nominal exposure. Instead, 
we adopt a computationally much faster approach, which 
utilizes the refractive approximation. That approxima¬ 
tion can be justified at these large scales by the following 
argument. 

Consider a two-dimensional turbulent patch on the 
phase screen of characteristic size r. Radiation can in¬ 
teract with that patch in two related but distinct ways: 
First, if there is a gradient in the index of refraction over 
the patch, it will act like a prism, so that the light will 
be refracted at a finite angle. This displaces the image of 
a star, i.e. gives rise to image motion without increasing 
the image width. Second, the finite size of the turbulent 
patch gives rise to diffraction. This increases the width 
of the image, without, in general, causing significant dis¬ 
placement. 
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The two effects depend very differently on the size of 
the patch. As discussed earlier, the phase power for Kol¬ 
mogorov turbulence is proportional to 7*3. We show be¬ 
low that the image displacement due to refraction is given 

by 


so 69 is proportional to r _ s. For diffraction on the other 
hand, 69 is proportional to r^ 1 , so the former is dominant 
on large scales, while the latter is dominant on small 
scales. The two effects are nearly exactly equal at the 
scale of the Fried parameter, ro, which is a measure of 
the seeing, and is defined by the relationship 

r 0 = ^0.423fc 2 J dzC^(z) 

For good seeing sites, ~ 0.6 arcsec, r 0 ~ 10 cm or so. 
Therefore on scales significantly larger than this, approx¬ 
imating the effect of the turbulence purely refractively is 
appropriate. 

The argument for the refractive approximation is even 
stronger if the primary effects one is interested in involve 
PSF asymmetry and decenter. Since a given scale is sam¬ 
pled multiple times over the pupil plane during a finite 
length observation, the net effect of image displacement 
due to that scale is reduced by a factor of N~ 2 , where 
N is the number of samplings. Since N is inversely pro¬ 
portional to r 2 then N~^ is proportional to r, so the 
contribution to PSF asymmetry and decenter scales like 
r», i.e. the largest scales completely dominate. 

The arguments above suggest that the natural choice 
for K cr n is of order 27 t /ro- At smaller wavenumbers, re¬ 
fraction is dominant, and we can ignore diffraction com¬ 
pletely. At larger wavenumbers, diffraction is dominant, 
but the full effect of the turbulence can be modeled statis¬ 
tically. The fact that tq is much smaller than the aper¬ 
ture of the telescope for modern telescopes like LSST 
means that this approach works extremely well for our 
application. The scale of the Fried parameter is sampled 
many times in the pupil plane, so even invoking the re¬ 
fractive approximation will merely cause circularly sym¬ 
metric image blur with a characteristic width of A/ro- 
That is what a true diffraction calculation would give 
anyway for the contribution from this scale. 

The analysis above suggests an especially simple imple¬ 
mentation approach for an image simulation code. The 
refractive approximation, formally, corresponds to an ex¬ 
pansion of as a Taylor series in r, and only keeping 
the first term 



4>{r) ss 4>o + f • V</> + ... 

When that substitution is made, we get 

PSF(n) = T(n — v—V</>) 

Z7T 

where T(n) is the diffraction profile of the telescope it¬ 
self. So the effect is purely image displacement. That 
means that the refractive approximation can be trivially 


implemented using ray optics: We bring photons individ¬ 
ually down through atmospheric layers, and simply de¬ 
flect them according to the local gradient of the phase. 
This will work as long as the resolution of the phase 
screens does not exceed K cr j^. 

This functionally accounts for the contribution 
PSF < (fh) As shown above, the total PSF is the convolu¬ 
tion of the PSFs from the lower and higher wavenumber 
regimes. The convolution of two probability distribu¬ 
tions gives the total probability of displacement if the 
contribution of each process is completely independent. 
For PSF > (rh ), we have a closed form expression for this 
probability distribution. We can thus sample that dis¬ 
tribution, in both magnitude and direction, and give the 
photon a second deflection accordingly. This two-kick 
raytrace then fully accounts for the net affect of that 
phase screen on all spatial scales. 

Numerically, we implement the two-kick raytrace by 
first generating two sets of turbulent screens (4 for each 
layer): the first (</><(r)) with turbulence generated from 
the van Karman spectrum only with wavenumber below 
Kent and the second (</>>(r)) with wavenumbers above 
K cr it • In practice, we found choosing K crit = 2n/1.5 r 0 
to produce the most stable and accurate results, which 
is consistent with the estimate discussed above. For the 
interactions of photons with the first set of screens we 
use the refractive approximation implemented by taking 
the derivative of the phase. This interaction produces 
the ellipticity and image motion of the PSF. The inter¬ 
action with the second set of screens is accomplished by 
asymptotically averaging the aperture-weighted Fourier 
transform of these screens and producing a single time- 
averaged PSF. This is equivalent to the long-exposure 
PSF for just this high frequency set of screens, and is 
equivalent to the closed-form solution we discussed above 
but in practice we use our exact generated screens. The 
effect of this second set of interaction adds additional 
wings to the PSF as expected. The two are implemented 
in sequence by performing the convolution in a Monte 
Carlo. We scale the angular displacement in propor¬ 
tion to A - 5 , which is consistent with the overall scal - 
ing from complete diffraction calculations (Fried 1965). 
If we choose n cr it to be an arbitrarily large number and 
therefore use the refractive approximation with the entire 
spectrum, then we get an unphysical circularl y symmet- 
ric PSF which is consistent with behaviour in deVries etl 


al. (2007), so the hybrid approach is necessary for proper 
treatment of the high frequency power. 

Thus in this approach we are making two approxima¬ 
tions that are testable: 1) that the refractive approxima¬ 
tion is valid for the first set of screens and any speckle 
inference structure averages sufficiently and 2) that the 
time-averaging of the small spatial scales is applicable. 
We test this combination of these two approximations 
quantitatively in §4 demonstrate its validity for reason¬ 
able exposure times t > Is. The net effect of this hybrid 
approach is that the computational efficiency is about 
four orders of magnitude larger than by using a pure 
Fourier approach (j 


5 A full diffraction calculation with no approximations would 
involve a series of mathematical operations over every screen pixel 
having size, p <C ro, for a screen area of Dvt exp where D is the tele¬ 
scope aperture and v is a typical wind velocity and t exp is the expo- 
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2.3.3. Photon Atmospheric Opacity and Dispersion 
Interactions 

We simulate molecular opacity by considering the 
wavelength-dependent optical depth through each seg¬ 
ment of the atmosphere layers. We use the HITRAN 
database to calculate the absorption line cross-section for 
the water and molecular oxygen absorption lines dRoth- 
man et al. 2009 Rothman et al. 2005 Rothman et al. 


20031 Kothman et al. [1998, Rothman et al.||1992[ IRoth- 


man et al.|1987 HITRAN|2014|). We follow the method¬ 
ology described in the appendix of Rothman et al. ( 2009) 
to calculate the absorption line intensity. We include 
pressure broadening so we get a dif ferent opacity a t each 
altitude. We use the formulae of IGreen et al.H988l for the 
calculation of the cross-section due to ozone. The opac¬ 
ity of the atmosphere is determined by calculating the 
local optical depth for each segment of the photons path. 
Between each atmospheric screen the column density of 
each molecular species of the atmosphere is calculated by 
integrating the density profiles in our bulk atmosphere 
model and the path length the photon. Thus, the prob¬ 
ability that the photon is destroyed along its particular 
path segment is equal to e~ T ^ where 

r = ^r l (A)=^ J <Ji(\,h)n,i(h)dh 


1 + 1CT 6 | 64.328 + 


29498.1 

146 


255.4 \ 

41 


Px 


1 + (1.049 - 0.0157T) 10" 6 P 6 0.0624 - 0 -°°° 2 680 

720.883 (1 + 0.003661T) (1 + 0.003661T)W 


where P is the air pressure (in mmHg), W is the water 
pressure (in mmHg), and T is the ground temperature (in 
C). This shifts the angle of the photon depending on its 
wavelength by a small angle. The direction of the shift is 
calculated by determining the vector to the zenith. The 
net atmospheric dispersion of the center of the field can 
be offset for a nominal wavelength. 


2.4. Telescope/Camera Optics Simulation 

The simulation of the telescope and camera optics is 
performed by defining a series of surfaces that separate 
the obvious volumes composed of various media (e.g. air, 
silicon, glass, vacuum). We first describe the geometry 
of the telescope and camera optics, and then discuss the 
photon interactions. 


2.4.1. Geometry of Surfaces and Media 

The optical elements of telescopes can often be de¬ 
scribed by cylindrically-symmetric aspheric surfaces, 


where the sum is taken over all the species of molecules. 
We also include the opacity due to aerosol scattering by 
including an additional contribution to the optical depth 
given by: 


T = T 0 


A 


500 nm 


where to = 0.02 and T = —1.28 are the default parame¬ 
ters, but vary at different sites due to a different mixture 
of aerosol types. 

When a photon hits a cloud screen pixel it has some 
probability of being absorbed (destroyed or lost by scat¬ 
tering outside the aperture). We use an average opacity 
of 0.85 magnitudes for each of the two clouds screens 
with a 1 a normal variation equal to the chosen mean 
opacity divided by 11.3. The latter value was chosen so 
that the structure function normalization is consistent 


with that determined by Ivezic et al. (2007). We chose 
to make the vari ance proporti o nal to ~the mean since the 
measurements of Ivezic et al. (2007) demonstrated that 
when the total cloud opacity is larger, there also seems 
to be a proportionally larger photometric variation. This 
results in an average total opacity of 1.2 magnitudes, and 
has a broad distribution, and is reasonably well matched 
to typical photometric variation of actual sites. These 
parameters may be different for different sites. 

Atmospheric dispersio n is simulate d by using the index 
of refraction in air from Filippenko (1982) 


sure time. Alternatively, our approximation only involves a math¬ 
ematical operation for each photon for the refractive calculation, 
and a one time diffraction calculation for all sources. Thus, the 
computational speed up scales as (DvtexpNsource)/(N p h 0 t 0n p 2 ). 
For D = 10m, v = 10m/s, t exp = 10s, N photon /N s ource = 1000, 
and p — 1cm, this ratio is 10 4 and is close to the actual observed 
computational ratio. 


z(r) = 


7 ? ^1 + — (1 + 


10 

n “if* 

i =3 


where the height as a function of radius, z(r), is expressed 
in terms of a radius of curvature, R, a conic constant, k , 
and higher order coefficients, a,;. Lenses are described by 
two aspheric surfaces (front and back), whereas mirrors 
are described by a single surface. The specification of 
these surfaces then define the regions where the glass 
is located. In those regions we comp ute the in dex o f 
refraction using the Sellmeier equation (Sellmeier 1871), 


n 



n 

A 2 -Ci 


B 2 A 2 b 3 a 2 

A 2 - C 2 A 2 - C 3 


We specify the detector plane by a series of rectangular 
devices each having a nominal center (x,y), a number (n x , 
n y ) of square pixels having fixed pixel size, p. For the 
index of refraction inside the Silicon, we use 


n = 3.36 + 0.211 E + 2.79e“tra5F- 


where E is the energy of the photon in eV (Phillip & Taft 
1960). Additionally, the regions of the telescope/camera 
system that are not h eld at vacu um ha ve the refraction 
index of air given by Filippenko (19821 as described in 
§2.3. We also model the three-dimensional spider sup¬ 
port structure as a series of rectangular volumes. There¬ 
fore, the combination of the series of aspheric surfaces, 
detector elements, and the indices of refraction of the 
media uniquely specify the configuration of the tele¬ 
scope/camera optics design. 

An important aspect of the simulation is to perturb 
the optical elements and detector segments from their 
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nominal positions due to expected small misalignments. 
With real telescopes a large part of the image quality 
budget is dominated by errors from misalignments and 
surface deformations of the optical elements. This is due 
to either fabrication or assembly errors, or environmental 
effects. Environmental effects include changes in temper¬ 
ature that misalign or deform optical elements, gravita¬ 
tional vector changes due to the different pointing ori¬ 
entations inducing misalignments and deformations, and 
other mechanical forces including wind and seismic vari¬ 
ations. Modern telescopes use control systems that at¬ 
tempt to correct for most of these environmental changes, 
but also may induce additional misalignments and sur¬ 
face deformations due to their residual inaccuracy. To 
account for all of these possible perturbations, the sim¬ 
ulator includes a misalignment and surface deformation 
for every optical element. The misalignments include the 
6 degrees of freedom: decenter (2), defocus, and Euler ro¬ 
tations (3). We model the surface perturbations using a 
Zernike expansion up to 5th order. Parameters describ¬ 
ing the tolerance for each degree of freedom are input to 
PhoSim for any optical design. We are currently devel¬ 
oping an interface to input actual engineering finite ele¬ 
ment calculations for the distortions and misalignments 
of surfaces for the LSST system specifically. The simula¬ 
tor can then choose a distorted surface based on various 
state variables (temperature, elevation, and actuator po¬ 
sitions) . 

A tracking model perturbs the entire telescope and 
camera system. The model simply perturbs the posi¬ 
tion of the photons in the reference frame of the camera 
and telescope, and represents the residual tracking errors 
that are expected for a nominal tracking system with 
parameters describing the tracking tolerance and update 
timescale (default of 0.1 seconds). We have implemented 
a Gaussian random walk model that varies throughout 
the exposure sequence. Every 0.1 seconds, a random 
walk step is taken in both elevation, azimuth, and rota¬ 
tion of the camera. The mean step is calculated so that 
the final RMS size of the jitter after a given exposure 
time meets the expected tolerance. Thus, the temporal 
spectrum is purely white up to 0.1 seconds. Between ev¬ 
ery 0.1 seconds, the jitter is interpolated. We also include 
errors in the effective exposure time of the shutter. 

2.4.2. Photon Interactions with the Optics 

As the photon is propagated through the tele¬ 
scope/camera optical system, it experiences a set of pho¬ 
ton interactions. The first involves dome seeing. We 
model dome seeing by perturbing the photon’s trajec¬ 
tory by an isotropic Gaussian angular distribution equal 
to the expected contribution. If the effective eddy size in 
the dome is sufficiently small this is a valid approxima¬ 
tion, but if it is not then some of the dome turbulence 
can still be accommodated in the ground layer. Some 
of the dome turbulence may have a limited drift speed 
depending on the airflow in the dome. 

An essential calculation is to find the location of each 
photon hit on a given surface. Consider a photon at 
position (x, y, z) with a unit vector trajectory of h. We 
loop through all possible surfaces by calculating the ray 
intersection distance. In order to find the intersection of 
a ray with a given surface with height, 2 = f(x,y ), we 
move the ray a scalar distance, l , and minimize, <5 


(5 =| z + h z l — z{x + h x l, y + n y l)\ 

where the surface includes the Zernike deformations. Be¬ 
fore calculating the intercept we make a three dimen¬ 
sional Euler transformation and spatial transformation 
of the photon’s position and trajectory to place the pho¬ 
ton into the frame of the optic. The above equation can 
be solved exactly (<5 = 0) for quadratic surfaces, so we 
first approximate z(x, y) with a parabolic approximation, 
z(x,y ), and solve the equation for l = lo. We then com¬ 
pute 6 between the actual surface and the position using 
the propagation distance, Zo, and choose a new value of 
l equal to 


* — HJ V * . 

n z 

We iteratively change l using this equation until it con¬ 
verges to within a tolerance of 0.01 microns. For even 
highly aspheric surfaces the method usually converges in 
3 to 5 iterations, which is essential for the computational 
efficiency of the code, since there are many surfaces in a 
typical telescope. 

When the next interaction surface is chosen, the pho¬ 
ton is moved to the new position. If the surface is a 
mirror we reflect the ray across the normal to the mir¬ 
ror using a three-vector computation. The normal has 
been pre-calculated for all surfaces including the surface 
perturbations. The ray’s propagation vector, h is then 
modified. If the surface is a lens, we refract the ray by 
applying Snell’s law, also using a three-vector compu¬ 
tation with the indices of refraction on the two sides of 
the surfaces determined by the formulae discussed above. 
We currently assume that the index of refraction is con¬ 
stant within the material for a given wavelength, but 
implementation of spatial varying indices of refraction 
is straight-forward. The detector elements and lens el¬ 
ements are treated identically as equivalent optical sur¬ 
faces. 

The optical elements have an interference coating that 
may affect the transmission of the photon. The proba¬ 
bility of transmission is a function of both incident an¬ 
gle and the wavelength. We include the full wavelength 
band in the simulation, so out of band filter leaks can 
be properly modeled. The two dimensional transmission 
probability is calculated using a full electro-magnetic in¬ 
terference boundary calculation through the actual sur¬ 
faces. When the photon reaches a coating, we use its 
wavelength and angle to decide whether it is reflected or 
transmitted. If the interacting surface is the mirror, the 
photon is destroyed if it is not transmitted. In the case 
of filter coatings or lens and detector anti-reflective coat¬ 
ings, the photon can be reflected backwards. Thus, we 
use the same reflection algorithm and allow the photon 
to propagate backwards. This implies that ghost pat¬ 
terns are included in the simulation. We currently have 
a four column interface which accepts coating reflection 
and transmission functions that are a function of both 
angle and wavelength. Coating descriptions are typically 
defined in this form, but we also have an external code 
for the complete EM multi-layer calculation which can 
calculate this format when the multi-layer structure is 
known. 



12 


Peterson et al. 


Note that we have already included the effect of the 
diffraction of the telescope pupil in the small-scale phase 
screen of §2.3.2. The point-spread-function therefore 
properly includes an Airy-like component due to the en¬ 
trance pupil. There is some freedom, however, in whether 
or not to imprint the spider pattern in the Fourier calcu¬ 
lation. If we do that, there is a significant variation in the 
projected spider size when light at large off-axis angles 
is included and, in principle, every source has a slightly 
different diffraction pattern. Alternatively, we can use 
the edge diffracti on calculation method of |Freniere, Gre-| 


gory, & Hassler (1999), where the photon’s position is 
shifted by an angular deflection of A/ (47rd), where A is 
the wavelength of the photon and d is the closest dis¬ 
tance a photon ray gets to the edge of any part of the 
spider structure. Thus, d can be calculated in fully three- 
dimensions and this calculation then results in both the 
correct geometric shadowing of the spider structure as 
well as the radial envelope of the diffraction spikes, but 
not any interference modulation of the diffraction spike 
pattern. 

To represent the incoherent scattering that occurs from 
micro-roughness on the mirror surfaces, we use a simple 
empirical model for large angle scattering. The micro¬ 
roughness of mirrors (at the nm level) primarily causes 
photons to scatter to very large angles (few arcminutes). 
At the current time, we have not implemented a physi¬ 
cal model for this, but instead invoked an empirical ra¬ 
dial distribution determined for stars measured with the 
Gemini South telescope: 


1 



We set the fraction of light in this diffuse halo compared 
to the core at a fixed fraction, /. Therefore, at the start 
of the telescope simulation the photon has a probability, 
/, of being scattered according to the above formula. The 
best fit profile had a value of (/, n, fo) = (0.135, 3.5,0.1°). 
There may be a contribution due to Mie scat terin g from 


dust particles included in th e above formulae (|King 
Roddier|[l995 Racine|[l996 ). 


1971 


To determine the probability of photoelectron conver¬ 
sion in the silicon detectors, we calculate the mean free 
path of the photon as it enters the silicon. The absorp¬ 
tion coefficient depends on the device temperature and 


the ph oton energy, acco rding to the model of Ra- 
jkanan, Singh, & Shewchun (1979), 


a(T) 


J2 C * A * 

i=l,2]j=l,2 


(E 1 — E g j(T) + E pi ) 
ew — 1 


(E 1 E g j(T) E p i) 

1 — e _ vr 

where E g {T) is given by 

BT 2 

E g (T) = E g ( 0)-^ 

where /3 = 7.021 x 10~ 4 eVK-\ 7 = 1108R, E gl (0) = 


+ A d \J E 1 — E g d(T) 


1.1557eV, E g 2 (0) = 2.5eF, E gd ( 0) = 3.2eF, E pl = 
1.827 x 10” 2 eF, E p2 = 5.773 x 10~ 2 eV, C x = 5.5, 
C 2 = 4.0, Ai = 3.231 x 10 2 cm- 1 eR- 2 , A 2 = 7.237 x 
10 3 cm _1 eV~ 2 , and A d = 1.052 x 10 6 cm _1 eV~ 2 . The 
mean free path is then calculated by taking the inverse 
of the absorption coefficient. We then calculate the ac¬ 
tual conversion path length by taking the mean free path 
and multiplying by an exponentially distributed random 
number. 

If the conversion length exceeds the full depth of the 
silicon, we allow for reflection off the backside of the de¬ 
vice. If reflection occurs, we continue the conversion cal¬ 
culation on the reflected ray. The reflection probability 
actually depends on a full EM interference calculation 
that can result in fringing. We calculate the reflection 
probability using a single layer of silicon with a height 
that is a function of position that depends on our per¬ 
turbations using an EM single-layer calculation. The sil¬ 
icon interference probability then depends on the index 
of refraction on the front surface and on the back surface 
as well as the photon polarization. 

We also include the possibility that the photon con¬ 
verts in an effective field-free region at the back surface. 
On the back surface of devices there may be a resid¬ 
ual field-free region due to the manufacturing process 
(such as using laser annealing). To simulate this effect, 
we simply remove any electrons that have converted in 
the field-free dead layer. This then will have the correct 
wavelength-dependence based on the photon conversion 
mean free path. 

2.4.3. Electron Interactions 

After the photon has converted to an electron, we sim¬ 
ulate the charge diffusion profile as it is drifted to the 
readout. To do this, we have developed a model of the 
electric field profile in the silicon. 

v a r 

E z {z) = -- 1 -/ dzrid(z) 

tSi Jtsi 

where V is the overdepletion potential, tsi is the silicon 
thickness, est is the permittivity in silicon, and n d is the 
doping density function which is given by 


_^ 

n d {z) = n bu i k + n b e a b +rife s f 

Note, that the impurity density, ribuik is not necessarily 
a constant due to the difference in segregation coefficient 
between the dopant and the silicon. The impurity may 
have a “tree-ring” pattern centered on the axis of the 
original boule. 

The relevant electron transverse diffusion at each 
h eight is calculated with Gaussian diffusion width, 
Dt c , where D is the diffusion coefficient, D = 
(fj, q (E,T)kT) /q, and the collection time is t c = 
f~ | g [ et)e (-)l ' We ^ ave deluded the effect of the ve¬ 
locity saturation of the electron in the expression for 
dg(E,T) 


A tq{E, T) — 


1.53 x io 9 T-°- 87 


L.01T 1 - 55 (l + E/ (l.OlT 1 - 55 )' 3 ) 
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where /? = 2.57 x 10~ 2 T°- 66 , T is given in K, and E is 
given in V/cm. After we compute the charge diffusion 
at the position in the silicon, we use the Gaussian width 
to move the electron laterally in the silicon and place it 
at the readout surface. Finally, the electron’s position is 
quantized by determining in which pixel it is located. 

In addition to the diffusion there is a small lateral mean 
shift due to any small lateral field. These lateral field re¬ 
sult from impurity variations, edge effects, charge stops, 
and accumulated charges during the exposure. Charac¬ 
terizing these lateral fields_in sili con devices is an active 
area of research (e.g. Kotov et al.|2006). The lateral kick 
the photon receives during its drift is given by 


Aa; = 

j dz 

Ay = 

1 dz 


E x {x,y,z) 

E z {x, y, z) 

E v {x,y,z) 

E z(x,y,z) 

We allow for simple parameterizations of the transverse 
field to add this complication. We are continuing to eval¬ 
uate the ideal parameterization based on real devices. 
The lateral field is also known to mod ified by the accu¬ 
mulated charge (Antilogus et al. 2014). In addition to 
lateral field, real devices may not have perfectly square 
pixels due to lithography errors, which can be simulated 
by simply making a non-regular map for the pixels at the 
readout. 

The simulation proceeds by collecting electrons in pix¬ 
els. We simulate the effect of charge saturation and 
bleeding by first not allowing a given pixel to exceed the 
full well depth in electrons, w. Once the full well depth 
is exceeded we move that electron towards the end of the 
row in either direction. We do not allow the electrons 
to move past the implant at the center of the device, 
and place the electron in the closest unfilled pixel along 
that row. Once the entire row has exceeded the full well 
depth, we remove the electron entirely. This then ap¬ 
proximates the effect of bleeding. 

Actual images of cosmic rays are added to the simu¬ 
lated im ages using real data from thick silicon devices 
Doty (priv. comm.). To do this, we constructed postage 
stamp images of 130 different actual cosmic ray events. 
We then add these randomly to our simulated images us¬ 
ing two important calculations. To determine how often 
to place a cosmic ray in the image, we use the produc¬ 
tion rate of 0.04 cosmic rays per sq. centimeter of silicon 
per second. The actual cosmic rays are expected to be 
a combination of gamma rays from local ground radi¬ 
ation and muons and other particles from atmospheric 
particle interactions. Our data have some combination 
of the two, but perhaps not in the correct proportions. 
A second calculation gives us the scaling of electrons in 
the cosmic ray data, to the appropriate volume of sili¬ 
con of the simulated pixels. This correctly normalizes 
the correct number of ionized electrons in the simulated 
devices. 


2.5. Electronic Readout Simulation 

Hot pixels are added to the image by randomly choos¬ 
ing a fraction of the pixels and then placing electrons 
equal to the full well depth. Similarly, a fraction of the 


pixels are flagged as dead and then the electrons are re¬ 
moved from those pixels. Hot columns are simulated by 
selecting some fraction of pixels that are the ends of hot 
columns and then setting those pixels and the pixels be¬ 
hind that pixel in the readout chain to full well depth. 
Dark current is computed for the length of exposure, and 
modeled by randomly adding a number of electrons to 
each pixel with Gaussian error. 

The CCDs are segmented according to the amplifier 
readout scheme. Rows or columns are added accord¬ 
ing to the number of pre-scan or over-scan pixels. Each 
pixel is then assigned a readout sequence according to 
the parallel and serial charge transfers. Electrons then 
have some probability of being shifted to a pixel behind 
it in either the serial or parallel direction during the read¬ 
out. We found that it was necessary to perform a shift 
for every pixel individually, since the CTE values can be 
quite high with modern devices (> 99.999%), and it is 
not possible to make a multinomial approximation for 
this effect. We then loop through the pixels in readout 
order. Read noise is implemented by using a Gaussian 
with mean equal to the expected read noise value. We 
also vary this value between amplifiers. 

Finally, the digitization process can be approximated 
by 


ADU = 


e 


+ B 


where e is the number of electrons in a given pixel, G is 
the gain, W is the full well depth, B is the bias, and IV is a 
non-linearity factor. We have not implemented a detailed 
model of gain and bias variations across each segment or 
as a function of time, but we do vary these parameters 
between amplifiers. We also have implemented possible 
digitization errors in which each bit is modified by 


biti = 


(ADU 

V 2‘ 



mod 2 


and then the final ADU value is given by 


ADU = biti2 l 

i 

where the cq values are empirically determined estimates 
that keep the digitization from being perfect. 

3. SOFTWARE IMPLEMENTATION 

The software is constructed to simulate a series of im¬ 
ages in identical form to what a real telescope would 
produce. The input resembles the combination of opera¬ 
tional commands a real telescope operator would execute, 
and a set of astrophysical catalogs. The final output is a 
sequence of FITS images produced from individual am¬ 
plifier chains of CCD devices. Following, we discuss the 
overall architecture and numerical implementation of the 
physics we described in §2. 


3.1. Architecture 

PhoSirn is written in object-oriented C++ code. The 
codes are run with Python scripting. The overall archi¬ 
tecture is shown in Figure 1. The C++ code is divided 
into 5 parts: the atmosphere creator, the instrument con¬ 
figuration, the trim program, the photon raytrace, and 
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the electron to ADC codes. The codes are configured 
to simulate a particular visit (a series of exposures at a 
given place on the sky). The first two codes set up all 
the input data that are required to describe the current 
atmosphere and instrument configuration. The trim pro¬ 
gram is then run for every chip where the astrophysical 
catalog is reduced to only objects that have a significant 
chance of producing photons on that particular chip (ei¬ 
ther object centered in the projected sky tile for that chip 
or particular bright objects that may have a large scat¬ 
tering halo or scattered ghost photons). This facilitates 
parallelizing the calculation. The photon raytrace sim¬ 
ulates the individual photons through the atmosphere, 
telescope, and camera and collects the converted photo¬ 
electrons in an image. Finally, the electron readout is 
simulated and the image is digitized in the electron to 
ADC code. 

The input and output to the code are well defined. 
The input consists of an instance catalog, which is a list 
of objects in the sky at the particular time of the obser¬ 
vation and a description of all the properties needed in 
§2.2. The instance catalog also includes the commands 
that a telescope operator would have available and other 
environmental parameters that may affect how the obser¬ 
vation would be done. Other astrophysical information, 
such as the position of the Sun and Moon is included as 
well. We also have a physics command file as an optional 
second input. This includes any commands to override 
our representation of the most realistic physics. This can 
be used in a large number of ways including by turning 
off a subset of effects in a modular way or setting parame¬ 
ters to specific values. Those are useful options, both for 
validation and testing, and also for studying the physics 
that might lead to a systematic error in a particular im¬ 
age processing algorithm. 

The main output of the code is raw digitized FITS 
images for every amplifier on every CCD. There are al¬ 
ternative other outputs as well that are unavailable with 
a real telescope. We can output an event file, which de¬ 
scribes the interaction position of every photon as it is 
propagated through each layer or surface. We also can 
output the actual number of photons detected from every 
source and the mean coordinate of those photons. This 
information would be only approximately known from 
the images. We also can output the relative throughput 
of photons at each layer or optical surface. 

The PhoSim code only has t wo package dependencies : 
cfitsio (Pence||l999) and fftw3 (|Frigo & Johnsonj[2005|), 
and is otherwise built using both standard C/C++ li- 
braries and custom numerical codes. This allows us 
greater custom numerical detail, and makes the installa¬ 
tion straightforward. The entire phosim code is designed 
so it can be implemented easily on grid computing. The 
architecture was designed so that I/O would be minimal 
and the simulations could be done in parallel at the chip 
level. The package can be run with both a script for 
laptop/desktop simulations and another script that uses 
CONDOR to generically r un s imulations on grid com¬ 
puting systems (CONDOR 2015). 

We have built an extensive validation framework with 
the PhoSim code. The framework includes both unit 
testing and integration testing. The unit testing executes 
individual functions to assess whether the return values 
obtain the correct values. Integration tests run the entire 


suite of code and determine whether measured properties 
of images obtain measured values within a specified toler¬ 
ance. Integration tests use instance catalogs of a limited 
number of objects (usually stars and galaxies), and then 
use the configuration files to run the photon simulator 
in a variety of configurations. We describe some of the 
results of the validation tests in §4. 

The entire Photon Simulation package is 

on a Git repository that is available at 

https: //www.bitbucket.org/phosim/phosim_release. 
Associated documentation is available at this site. We 
release a tagged version several times per year. We 
use a modern modular object-oriented software design 
approach where the speed of the code is a very important 
consideration that we discuss in the following section. 

3.2. Numerics and Optimizations 

There are a variety of numerical implementation de¬ 
tails specific to each of the implemented physical effects. 
In general, Monte Carlo simulation times are propor¬ 
tional to the number of points in the Monte Carlo in¬ 
tegration (in our case, the numbers of photons) and a 
fixed overhead associated with setup. Either reducing 
the time per photon or reducing the number photons that 
need to be simulated can minimize the simulation time. 
For the former, minimizing the number of mathemati¬ 
cal operations done on each photon reduces simulation 
time. Removing redundant calculations by saving values 
in pre-calculated tables wherever possible is one key to 
optimization. We do this for the shapes of the optical 
surfaces, transmission curves, turbulence screens, optical 
depths, etc. The overall reduction of the number of lines 
of code that need to be executed in the inner loops is 
also important, so this has been a priority throughout 
its development. Currently, we are simulating a photon 
in about 2 ps on a 2.5 GHz processor, implying a very ef¬ 
ficient number of calculations per photon. This can pos¬ 
sibly be improved further (and has during the recent de¬ 
velopment), but it is a fairly small number of operations 
considering the detailed physics and fidelity constraints. 
We profile the code periodically, and the various parts of 
the calculation described in §2 contribute roughly equally 
at the current time since obvious bottlenecks have been 
removed. 

We have performed another set of optimizations that 
reduce the overall numbers of photons that need to be 
simulated. These optimizations have a minimal effect on 
the fidelity, and all the optimization can be turned off 
for detailed comparisons. We have three such optimiza¬ 
tions that we call: 1) dynamic transmission optimization, 
2) saturated object simulation, 3) and background op¬ 
timization simulation. Dynamic transmission optimiza¬ 
tion works by attempting to guess whether the given pho¬ 
ton has any chance of surviving the Monte Carlo simula¬ 
tion. The most common way for a photon to be removed 
is by not surviving the filter transmission. A simple op¬ 
timization is to pre-roll the random numbers at the start 
of the simulation of the photon. We then store a worst- 
case transmission curve for each atmospheric layer and 
surface coatings for each wavelength as the simulation is 
running. This is necessary because of the complication 
that our transmission functions are not only wavelength- 
dependent but they may be angle-dependent and time- 
dependent as well. Therefore, we simply estimate first if 
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the photon has any chance of surviving the series of trans¬ 
mission functions before the photon is simulated through 
the full physics. This improves the simulation speed by 
an order of magnitude, and has a negligible change on 
the photometry. 

A second optimization is the saturated object simula¬ 
tion optimization. The brightest stars in a typical field 
contain a very large fraction of all the photons. However, 
the vast majority of those photons end up in a saturated 
bleed trail that has a negligible amount of useful infor¬ 
mation. To take advantage of this, we have built an opti¬ 
mization that makes the ray represent N photons instead 
of just one when it is highly likely that ray will end up 
in a saturated bleed trail. This can be done at run-time 
as soon as the pixels start to saturate and exceed the full 
well depth for a given source. To simulate the wings of 
the saturated source accurately with the right statisti¬ 
cal properties, we are able to keep track of the rays that 
might end up in the wing from large angle scattering 
(either diffraction spikes, mirror surface microroughness, 
or diffraction in the atmosphere). We can then enhance 
the probability of those events happening by artificially 
looping over that physics N times. We then choose N 
such that the probability of the photon being kicked by 
an angular distance greater than the current minimum 
radius of the saturation pattern, ro, is greater than a 
value, a. Normally, a would be a small quantity (few 
percent) without optimization, but here we enhance its 
probability to 90%. We preserve probabilities, however, 
by letting the ray represent M photons, if it is not a large 
angle photon. M is given by 


Table 1 

Required Instrument and Site Specific Data 


Type of Data 

Information 

Optical prescription for 
every optical element 

type of optic, R CU rv 
position, outer/inner radii, k, 

£* 3 , £* 4 , 015 , £*6 , a 7, <* 8 ? <^ 9 , <^10 
medium adjacent to surface 

Focal plane geometry 
for every device 

position, pixel size, pixel #’s, 
thickness, readout speed, 
shape distortion parameters 

Surface perturbation 
and misalignments for 
every degree of freedom 

distribution parameters 

Coating for 
every optical surface 

transmission probability 
vs. angle and wavelength 

Device data 
for every sensor 

hot/dead pixel rate, 
hot column rate 


Readout Data amplifier segmentation, 

for every amplifier pre/over-scans, read noise, 


serial and parallel CTE, gain, 
bias, non-linearity 


Site data 

turbulence distribution 

for every layer 

outer scale distribution 
wind vector distribution 

Location 

latitude, longitude, altitude 


not simulate many more than ^fp where p is the pho¬ 
tons per pixel. This results in two orders of magnitude 
faster background simulation, and varying degrees of ac¬ 
curacy depending on the choice of N and er. With the 
default parameters, a reduced y 2 of 1.1 is obtained for 
comparing the images of off-axis chip simulation with this 
optimization turned on and off. For detailed studies of 
either bright stars or background models the optimiza¬ 
tions can be turned off so photons are properly simulated 
one photon at a time for complete accuracy. 


On the other hand, we only let the ray represent one 
photon if it is a large angle photon. Similarly, we con¬ 
serve the photon detection probability by letting the ray 
represent N photons, if it is removed. Thus, the algo¬ 
rithm still simulates one photon at a time for virtually 
all photons that will be measurable in the image (the 
wings), but simulates several at a time for photons in 
the saturated bleed trail or those not detected at all. 

The final optimization involves the simulation of the 
background photons from airglow, scattered moonlight, 
twilight, or dome light simulations. Since these photons 
outnumber the astrophysical photons, it is important to 
simulate them efficiently. On the other hand, we found 
that using parametric models of the background (i.e. try¬ 
ing to predict the flux in each pixel and then adding 
noise) did not result in high enough accuracy for some 
of the more subtle physics details. An algorithm that 
facilitates faster simulations is to represent the ray as 
N photons for part of the simulation. When the pho¬ 
ton is closer to the pupil plane, a diffuse illumination 
pattern will produce an nearly identical photometric re¬ 
sponse for small angular distances. Equivalently, most of 
the differences in background from one pixel to the next 
occur because of physical effects near the image plane. 
Therefore, when the photon reaches the detector we can 
randomly spread the N photons in a Gaussian pattern 
(with a of several arcseconds width) and simulate the 
detector physics one photon at a time. To not induce 
fluctuation patterns, N should be chosen so that we do 


4. RESULTS 

We implemented the LSST design details through a se¬ 
ries of data input files. The physics code is written delib¬ 
erately without including any reference to LSST-specific 
data, so implementation of other telescopes is straight¬ 
forward through an alternative set of design files. We 
do not describe the detailed design parameters of LSST 
here, but Table 1 lists the variety of design information 
that is needed to describe a telescope and site. The 
physics implemented above should be appropriate for 
most optical survey telescopes without significant modi¬ 
fication. There may be possible extensions to other spe¬ 
cialized systems as well, such as an adaptive optics tele¬ 
scope, where the effect of the AO control loop on the 
residual phase error would have to be considered. An 
application to a space-based telescope would be straight¬ 
forward as well. Extending the physics wavelength cov¬ 
erage into the UV and further into the IR would also 
require little modification. However, the validation stud¬ 
ies we pursue below are most appropriate for simulating 
optical survey telescopes. 

In Figure [5] we show the path of the photon through 
the large dynamic range of scales in the photon simula¬ 
tion. The top left image shows the path of photon in a 
cylindrical column through the atmosphere, the bottom 
left and bottom center images show the photons mov¬ 
ing through the telescope system, and the bottom right 
image shows the conversion of photo-electrons in the sil¬ 
icon. The top image then shows the resulting images of 
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Table 2 

Connection between relevant physics and image properties 


Image Property 

Most Relevant Physics 

Photometric Acceptance 

Geometric Acceptance, Coatings, 
Photo-electric Conversion, 
Atmosphere Opacity, Clouds 

PSF size 

Optical Design 

Atmospheric Turbulence 
Perturbations and Misalignments, 
Dome Seeing, Tracking 

Charge Diffusion 

Mirror Micro-roughness 

Astrometric Scale 

Optical Design, 

Atmospheric Dispersion, 

Tracking, 

Perturbations and Misalignments 

Background Level 
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accurately reproduce the image properties from a given 
telescope and site). Representational errors are much 
harder to address, since they involve understanding the 
exact characteristics of the telescope, camera, and site 
and not just the correct implementation of the physics 
of photon and electron propagation. 

The required calculational accuracy of PhoSim is a 
complex topic. As described in the previous sections, 
we avoid making approximations, unless they result in 
orders of magnitude faster simulation rates. Thus our 
goal is to make the physics as complete and accurate 
as possible. Although there is a vast range of science 
applications that would place disparate requirements on 
simulation accuracy, a first order estimate of required 
accuracy can be obtained by considering the statistical 
error for measuring various attributes of the point spread 
function. If the calculational error is significantly below 
the statistical error for a source with a given number of 
photons then the error will be unobservable. In the Ap¬ 
pendix, we estimate generically that in an optical survey 
the brightest new source will have statistical photomet¬ 
ric error of 10 millimags, a FWHM uncertainty of 12 
milliarcseconds, a centroid uncertainty of 7 milliarcsec- 
onds, and an ellipticity uncertainty of 1.4%. Below, we 
show that the known calculational errors of PhoSim are 
significantly below those thresholds. 

For the calculational validation of the atmosphere sim¬ 
ulation, we show typical phase screens in Figure[9] Three 
examples of 50m by 40m phase screens of the combined 
7 layer set of screens after combining the screens on the 
4 different pixel scales are shown. Despite the com¬ 
plexity in constructing these images, they qualitatively 


stars and galaxies are they are collected in pixels. This 
figure does not completely show the physics detail in the 
simulation. Figure [3j however, demonstrates the physics 
detail by simulating a single star and successively turn¬ 
ing on more physical effects in the simulation. Each sep¬ 
arate part of the simulation contributes in different ways 
to the size and shape of the PSF, the photometric in¬ 
tensity, and the astrometric position of the image. The 
star was simulated in the u,i,y filters and combined to 
be a three-color RGB image, so that the chromatic ef¬ 
fects can be seen. Figure [4] shows a simulation of a star 
through the same optics and atmosphere configuration, 
but at a point 1 degree away from the center of the field. 
Through comparison of these two figures, the subtle spa¬ 
tial dependence of the PSF and photometric properties 
can be seen. Table 2 summarizes the relevant physical 
effects that determine the particular image properties. 
Figure [5] is a collection of 3 amplifier images wi th various 
stars and galaxies u sing catalogs generated by |Connolly| 


jLane, Glindemann & Dainty 

1992 Ellerbroek 

Le Louarn 

2002 

Britton 

2004 

Jolissamt 2010), 


|et al.| (|priv. comm.). Every photon has been simulated 


and sampled from the spectral energy distributions and 
spatial models in the catalog. 

To date, we have validated the most critical aspects of 
the simulator. However, a complete validation is beyond 
the scope of this work and awaits wider community in¬ 
volvement using real data from astronomical telescopes. 
Following, we discuss the most critical aspects of address¬ 
ing any possible calculational errors (i.e. given an ideal 
setup can PhoSim properly calculate a set of given quan¬ 
tities). The larger future validation, however, involves 
possible representational errors (i.e. does the simulator 


2002 

but we have a representation o 
kilometers. A quantitative validation is shown in the up¬ 
per left panel of Figure [9] where we calculate the two 
dimensional structure function of the p hase screen and 
compare with the analytic calculation of Fried (1965) for 
a pure Kolmogorov spect rum and the numerical integral 
of Lucke & Young| (|2007D for various values of the outer 
scale for a von Karman spectrum. No visible artifacts 
can be seen in the structure function. The propagation 
of light through the phase screens uses the hybrid propa¬ 
gation technique that can be compared with a tra diti onal 
Fourier approach. This is shown in Figures [lO] and [XT] 
Fig ure 1 10 1 shows an instantaneous exposure, whereas Fig¬ 
ure |11| shows a typical 15 second LSST exposure. The 
PSF using both numerical approaches is calculated and 
compared on LSST pixel scale (bottom) and finer scale 
(top). The residual subtracted PSFs (right) show no ob¬ 
vious biases or errors. Quantitative comparisons can be 
made by measuring weighted moments and computing 
the centroid and ellipticity. Even for bright stars, the 
statistical errors in the centroid and ellipticity are statis¬ 
tically consistent using the two approaches. To quantify 
this, we performed two dozen simulations of stars using 
the two approaches and plot the measured ellipticities 
and centroids demonstrating a high correlation in Fig¬ 
ure |12| The differences are shown in the histogram in 
the right panel and are consistent within statistical er¬ 
rors. This histogram can be used to set an upper limit to 
the error of the technique of 0.60 milliarcseconds for the 
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centroid error and 0.0046 for the ellipticty error. Thus, 
the only significant error of this light propagation tech¬ 
nique may be for arbitrarily short exposures with diffrac¬ 
tion limited telescopes where the h ybri d technique will 
not reproduce the speckles in Figure [lOj However, as we 
discussed in §1.3, this is not relevantror optical survey 
telescopes. 

For the calculational validation of the instrument sim¬ 
ulation, we show a spot diagram simulation in Figure [l3l 
A spot diagram is produced in PhoSim by turning off oe^ 
tector effects, ignoring perturbations and misalignments 
of the optics, and by making arbitrarily small pixels. We 
compare the result of 5 different monochromatic simula¬ 
tions for a 1.4 d egree of f-axi s point with the commercial 
raytrace code, Zemax (2014). The rms spot size is about 
0.1 arcsecond, so smaller than a single pixel. The x-axis 
in Figure [l3| subtracts off the large positional offset ex¬ 
pected for such an off-axis source of 254911.0 uni. This 
demonstrates that we are correctly predicting the posi¬ 
tions of photons at least to the 7th digit. A detailed ray- 
by-ray comparison can be used to assess the quantative 
accuracy of the geometric raytrace. We find that the dis¬ 
crepancy in ray positions has an average displacement of 
0.18/im, most likely due to the numerical accuracy with 
which we store surface maps. Note that this is less than 
l/50th of an LSST pixel. 

For the verification of photon throughput, we simu¬ 
late a top hat spectral energy distribution, and verify 
the analytic prediction for the number of photons. We 
then compare the detected number of photons with the 
analytic calculation 


1 m 500+48.6 o., (A 2 A]) 

photons = ylO ~ 2 - 5 texpTr (r 2 a - rj) -—- 

n r \ j 5 QQ nm 

where h is Planck’s constant, 771500 is the AB magnitude 
at 500 nm, texp is the exposure time, r 0 and are the 
inner and outer radii of the aperture, and A 2 and Ai are 
the wavelength limits of the top hat. Figure 14 shows 
the number of photons generated using PhoSim as com¬ 
pared with the analytic prediction using square spectral 
energy distributions, the aperture of LSST, and the ex¬ 
posure time. The results are consistent within statistical 
errors, and well below any reasonable science application 
with PhoSim since the mean inaccuracy is 0.45 millimags 
(0.045%). Thus, the numerical simulation of the atmo¬ 
sphere and instrument is sufficient for nearly all practical 
science cases. There is still considerable work to validate 
represent atio nal errors. 

Figure [15] shows the speed of the calculation, and 
demonstrates the efficiency of the photon Monte Carlo 
approach. The speed is proportional to the number 
of photons for unsaturated sources (above 15th magni¬ 
tude) and is approximately 400,000 photons/s on a typ¬ 
ical workstation (2.5 Ghz Mac Intel Processor). Back¬ 
ground photon simulations are considerably faster (fac¬ 
tor of 60), and saturated sources are as well (factor of 8 
for a 12th magnitude star) due to the optimization de¬ 
scribed in §3.2. The computation takes about 3 hours 
for a single chip of LSST (13 by 13 arcminutes) filled 
with a typical distribution of stars and galaxies for a 15 
second exposure. Thus, despite the physics complexity 
pursued in this work, we can simulate typical individ¬ 
ual sources which may have thousands of photons in a 


few milliseconds. Using large scale computing, we have 
run the code on about 2000 processors simultaneously 
for a variety of data challenges to simulate LSST simula¬ 
tions for the data management team. Therefore, during 
these runs we are already generating data within reach of 
the actual real time image production rates of the future 
highest data rate telescope (LSST) j^] Since the simula¬ 
tions are done in parallel, there is no barrier to scaling 
to arbitrarily large numbers of cores. 

5. CONCLUSION AND FUTURE WORK 

We have demonstrated the implementation of atmo¬ 
spheric and instrumental physical effects appropriate to 
astronomical image simulation for optical survey tele¬ 
scopes. The simulation approach using a photon Monte 
Carlo is both efficient and flexible. Simulations with this 
level of physics detail are remarkably fast with the pho¬ 
ton Monte Carlo methodology. Detailed second order 
image properties can be studied in analyses that require 
high precision (e.g. astrometry, detailed PSF modeling, 
photometric calibration) using this simulation tool. 

There are a variety of ways to improve and extend the 
code further. Future work involves validating physics 
techniques and adding more details to the representa¬ 
tion of the instrument and site characteristics. The most 
important fidelity improvements are: 1) a physical model 
predicting the perturbations and misalignments, 2) a de¬ 
scription of lateral charge diffusion associated with elec¬ 
tric field distortions, and 3) more accurate atmosphere 
site models for both turbulence and opacity. Detailed 
validation may be achieved by greater community in¬ 
volvement with the simulations of existing telescopes in 
various configurations, an we expect that this tool will 
be useful for many different science applications. 
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APPENDIX 

STATISTICAL ERROR OF AN OBJECT’S ATTRIBUTES IN AN OPTICAL SURVEY 

To determine the statistical error on a source’s attributes, we simulated a ellipsoidal gaussian with N photons with 
PSF of size, a. We then measured the source’s attributes (FWHM, ellipticity, flux, and centroid) by a weighted second 
moment method. We confirmed the scalings that the statistical error for the photometric flux is ./V -0 ' 5 , the error 
on the FWHM is 2 (21og2) 0 " 5 N~°- 5 a, the error on the centroid is 2°- 5 IV _0 ' 5 a and the error on the PSF ellipticity is 

20.5^y—0.5 

Science applications can be done with objects of all brightnesses and therefore values of photons, N, but note that 
for an optical survey most of the interesting new sources will be near the detection threshold. This occurs when 

N F > s\/ Fpb 

where N is the number of photons from a source in a single frame, F is the number of co-added frames, p is the effective 
number of pixels the source is distributed over, b is the background rate per pixel per frame, and s is the source signal 
to noise. Generically, p is around 10 to sample the PSF but not reduce sensitivity, b is in the range from 10 to 10 3 to 
be sky-noise limited but not limit the dynamic range of a CCD up to the full well depth (typically 10 5 with modern 
devices), and a reasonable detection threshold would have s of 10. Therefore, N would typically be between 100F -0 5 
and lOOOF -0 - 5 for a source at the survey’s detection threshold. Conservatively, if we consider sources an order of 
magnitude brighter than the detection threshold, then the brightest new source in a survey would then have around 
10 4 F -0 - 5 counts in a single exposure or 10 4 T 0 ' 5 counts in the co-added exposures. For seeing of 0.5”, the statistical 
uncertainty of the object’s attributes in a single exposure would have a photometric error of 10 millimags, a FWHM 
of 12 milliarcseconds, a centroid error of 7 milliarcseconds, and an ellipticity error of 1.4%. Therefore, a reasonable 
set of requirements for the accuracy of an optical survey simulator is that any known calculational inaccuracy is 
below several millimags of photometric error, several milliarcseconds of FWHM and centroid error, and percent-level 
ellipticity error. This would be sufficient to not add a significant systematic simulation error on an object’s attribute 
for a source with less than 10 4 photons in a single exposure or 10 4 combined photons in a co-added exposure. In most 
cases, the simulation error would not occur in the same way in a series of simulated exposures so the single epoch 
requirement would be sufficient, but if it did, the error would, for example, result in a factor of 10 stricter requirements 
for 100 frames. In either case, the scalings in this appendix can be used to estimate the errors on a given object that 
is the study of a science case. 
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Figure 1. The basic architecture of the phosim code. The five codes and the controlling scripts are shown as green hexagons. The inputs 
are the instance catalog, optional physics command files, and the static design data. The user can interact with the code running the PhoSim 
script and either construct instance catalogs, modify the physics commands, or change the input data. The integration validation tests 
are a combination of simple instance catalogs and physics commands. The main output is amplifier images, but a number of intermediate 
outputs are also available. 
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Figure 2. A graphical representation of the photon’s trajectories. The color of the photons is proportional to the photon wavelength. 
The upper left panel shows the path of photon’s through a cylindrical column of the atmosphere. The bottom left panel shows the 3 mirror 
LSST design. The bottom middle panel zooms in around the lenses and shows the non-blue photons being reflected backwards from the 
filter. The bottom right panel shows the photons (blue) converting into photo-electrons (black) as they drift to the readout. The top panel 
is a two dimensional histogram (image) of the final electron positions. 
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Figure 3. A simulation of a single star at the center of the field with different physics modules turned on in each of the 12 panels. The 
scale is different in each of the panels and is indicated by a bar representing 0.2 arcseconds. The simulation is done in three filters (u, i, 
y) through the same atmosphere and instrument configuration to show the chromatic effects of each physical effect. The PSF size/shape, 
photometric intensity, and astrometric position is therefore changed by each part of the simulation. The instrumental PSF is dominated 
by the perturbations, misalignments, and charge diffusion. The atmospheric PSF is produced by a combination of the ground layer and 
the free atmosphere seeing. 
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Figure 4. The same as the previous figure but with the star 1.4 degrees off-axis. Comparison of the two figures shows not only the 
complex nature of the simulation, but also the prediction of the complex field-dependence of these effects. 
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Figure 5. A simulation of the entire field (10 sq. degrees) of the LSST field of view. 10 million stars and galaxies are in the simulation 
with over 1 trillion photons. This simulation was executed using CONDOR grid computing for about 1000 CPU hours in which each 
individual CCD was simulated in parallel. The image has over 3 billion pixel, so the full detail cannot be observed. The variation in bias 
levels of the individual amplifier is visible, as well the vignetting of the background near the corners of the field. 
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Figure 6, 


The central 5x5 chips covering about 1 sq. degree. Individual bright stars are visible as well as hot columns. 
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Figure 7. A close up view of the central chip where more stars and brighter galaxies are visible. 




Figure 8. An image of 3 amplifiers (3/16th of a chip) showing the full view of the stars and galaxies that we have simulated. Every photon 
has had the full physical simulation described in this work. Typically, an average galaxy near the detection threshold (24th magnitude) 
only has a few thousand photons in a 15 second exposure. 
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Figure 9. Examples of three phase screens on 50m by 50m scales. The colors indicate the relative phase shift as a function of position. 
The images are generated by adding the phase screens from the 7 atmosphere layers and each phase screen on the 4 different pixel scales 
for each layer. No obvious artifacts resu lt from the numerical grids. The upper left panel shows compar ison wi th the structur e functions 
predicted for a pure Kolmogorov model < |Fried 1965[ ) and a von Karman model for different outer scales {Lucke &; Young|2007| > 
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Figure 10. A comparison of the point-spread-function induced by the atmosphere using a full Fourier approach (middle panels) and the 
two approximations discussed in the Appendix (left panels). The simulations are for a source exposed for 1.5 milliseconds. The top panels 
show the PSF before pixelization, and the right panels show the residual difference (colors represent the -3 (dark purple), -2 (purple), 
-1 (blue), 1 (yellow), 2 (orange), 3 (red) sigma residuals). The only visible differences are the speckle structure that occurs on l/20th 
of a pixel scale. The geometric approach PSF has FWHM of 0.800 =b 0.021, ellipticity of —0.113 zb 0.007, 0.082 =b 0.007 and centroid of 
0.015 =b 0.006, -0.003 =L 0.006. The fourier approach PSF has FWHM of 0.834 =b 0.021, ellipticity of -0.136 =b 0.007, 0.064 zb 0.007 and 
centroid of 0.013pra0.006, —0.004 zb 0.006. The reduced x 2 of the residual map is 1.36. 
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Figure 11. The same as the previous figure except for a 15 second exposure. The difference in the two simulation methods is purely 
statistical. The geometric approach PSF has FWHM of 0.0793 zb 0.021, ellipticity of —0.015 zb 0.007, —0.016 zb 0.007 and centroid of 
-0.007 =b 0.006, 0.003 =b 0.006. The fourier approach PSF has FWHM of 0.766 =b 0.021, ellipticity of -0.011 =b 0.007, 0.005 =b 0.007 and 
centroid of —0.009pm0.006, —0.001 zb 0.006. The reduced of the residual map is 0.84. 
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Figure 12. Measurements of the centroid and ellipticity of example atmospheres using the two techniques. The y-axis of the left plot uses 
the traditional Fourier light propagation technique, whereas the x-axis is the measurement using the geometric hybrid approach. There is 
a clear correlation demonstrating that both techniques are capturing the same details. The statistical differences are shown in the right 
panels. The distribution is consistent with expected statistical errors (red), and can be used to set an upper limit to the numerical technique 
(see text). 
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Raytrace PSF & Astrometry 
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Figure 13. Comparison of spot diagrams through the LSST design for an off-axis star at various wavelengths. Note the scale makes 
these images fill only 4x4 LSST pixels, so the detail is reproduced on a very fine scale. The x-axis has a large offset subtracted so this 
demonstrates that the rays are arriving at the focal plane to a very high accuracy. Detailed ray by ray comparison shows a rms difference 
of 0.018 microns, which is negligible for all practical science cases. 
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Throughput Accurocy Test 



Figure 14. The photon throughput using six different square spectral energy distributions in PhoSim vs. an exact analytic photometry 
integral (see text). The bottom plot shows the residual photometric error, and tests PhoSim’s ability to properly sample SEDs, accept or 
reject photons through various modules, simulate the correct geometric acceptance. 
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Speed Test 



Figure 15. The speed of the simulator. The diamonds show the number of photons from all astrophysical sources on a LSST chip (13 by 
13 arcminutes). The blue and red curves show how many stars and galaxies, respectively, per sq. degree. The yellow curve splits the sky 
background into equivalent sources placed every 15 arcseconds. The number of seconds for each magnitude bin is shown by the triangles. 
For faint non-background unsaturated sources the simulation time scales as the number of photons and is typically 400,000 photons per 
second. For saturated sources and background there are significant optimizations that increase the simulation time. A 12th magnitude 
source is simulated at about a factor of 5 faster than an unsaturated source, and the background is simulated about a factor of 50 faster 
than unsaturated sources. There is typically less than 1 source per LSST chip (13 by 13 arcminutes) below 12th magnitude, so those 
sources are not common. 













